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Abstract 

In  the  present  paper  an  implicit  time  accurate  approach  combined  with  multigrid  and  precon- 
ditioning is  used  for  the  large-eddy  simulation  of  low  Mach  number  flows.  It  will  be  shown  that 
the  present  approach  allows  an  efficiency  gain  of  a factor  4 to  7 compared  to  the  use  of  a purely 
explicit  approach.  The  efficiency  varies  according  to  the  test  case,  grid  clustering,  physical  time 
step  and  requested  residual  drop. 


1 Introduction 

Preconditioning  is  widely  used  as  a convergence  acceleration  tool  for  low  Mach  number  flows.  This 
technique,  however,  is  not,  time  accurate  excluding  its  straightforward  use  for  unsteady  problems. 
Multigrid  on  the  other  hand  is  one  of  the  most  efficient  convergence  acceleration  tools.  In  practice, 
most  multigrid  applications  are  for  steady  state  problems.  The  dual  time-stepping  approach  can 
be  used  for  time  accurate  problems  and  has  the  advantage  that  non-time  accurate  tools  such  as 
preconditioning,  multigrid,  residual  smoothing  and  local  time  stepping  can  be  used  in  a time  accurate 
context.  In  the  present  paper  the  use  of  these  methods  is  extended  towards  large-eddy  simulations 
of  low  Mach  number  flows. 

The  LES  code  is  an  extension  of  the  finite  volume  Reynolds  Averaged  Navier-Stokes  (RANS) 
solver  EURANUS,  [1].  EURANUS  contains  an  efficient  multigrid  solver  based  on  the  Full  Ap- 
proximation Storage  (FAS)  scheme  for  flows  ranging  from  low  subsonic  up  to  hypersonic,  [2]  [3].  The 
Runge-Kutta  scheme  acts  as  a smoother  for  the  multigrid.  For  low  Mach  number  flows,  the  efficiency 
of  multigrid  strongly  depends  on  the  stiffness  of  the  equation  and  the  ability  of  the  smoother  to  damp 
the  high  frequency  errors.  The  stiffness  is  removed  with  the  preconditioning  technique  which  was 
applied  successfully  to  different  steady  and  unsteady  RANS  flows,  [4] [2],  and  in  the  present  paper  will 
be  used  in  the  LES  context.  A central  second-order  finite  volume  spatial  discretization  is  used  in  the 
LES  code.  Fourth-order  artificial  dissipation  is  added  to  increase  the  high  frequency  damping  of  the 
smoother.  In  order  to  avoid  laminarization  of  the  flow,  the  artificial  dissipation  is  only  added  to  the 
continuity  equation  but  not  to  the  momentum  and  energy  equations.  According  to  our  experience, 
this  also  suffices  to  guarantee  a good  overall  convergence.  As  mentioned  above,  dual  time-stepping  is 
used  to  update  the  equations  in  time.  The  physical  time  derivative  is  discretized  either  by  a backward 
differencing  method  or  by  a trapezoidal  scheme.  Numerical  experiments  show  that  the  trapezoidal 
scheme  is  more  flexible  and  allows  bigger  time  steps  as  compared  to  the  backward  differencing. 

The  Smagorinsky  model  is  used  for  subgrid-scale  modeling.  Both  the  classical  model  and  the 
dynamic  model,  [7],  are  available  in  the  code. 
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2 Mathematical  Formulation 

The  preconditioning  procedure  in  the  framework  of  dual  time-stepping  can  be  formulated  as,  [2]: 
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Where  r is  the  pseudo-time,  t the  physical  time  and  <Ty,  T,j  respectively  the  filtered  stresses  and 
SGS  stresses: 


_ , dui  duj 

^ = /%-+£te- 


2 duk 
3dxkdij) 


Tij  = p(tl,Uj  - U,Uj  ) = -/J(( 

The  subgrid-scale  heat-flux  Qj  is  defined  as: 


dui  du j 

dxj  + dxi 


2 duk 

3 dxk 


(4) 

(5) 


Qi  — pCp(Tuj  Tv, ) (6) 

It  is  modeled  using  an  eddy  diffusivity  SGS  model: 


PtCp  dT 
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pt  = pCA2\S\,  and  Prt  = 0-5  or  is  calculated  dynamically,  0 and  a are  the  preconditioning 
parameters,  a is  taken  as  a constant  around  -1,  and  0 can  be  defined  locally  as,  [9], 


0 = MiniMaxi'/UJTi,  ■—),  c) 


(8) 


With  A the  smallest  cell  length,  l the  biggest  characteristic  dimension  of  the  domain,  At  the 
physical  time  step  and  c the  speed  of  sound.  In  the  present  LES  calculations  At  is  rather  small,  due 
to  physical  restrictions,  such  that  ^ is  the  dominant  term  in  (8).  As  a result,  0 will  be  constant 
on  the  whole  field. 

By  applying  the  finite  volume  method  to  (1),  (2)  and  (3),  they  can  be  written  in  matrix  form  as 
follows: 


Written  out  fully: 
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The  physical  time  derivative  is  discretized  with  a multi-step  scheme: 

r_jao$2  =_^LUn  + sn  + -riIies(u)  (11) 

or  A t 

whore  Sn  is  the  source  term  defined  as: 

5"  = ~(Unn)  ~ §T{f/n',0)  “ (»""*«)  +l2Res(Un)  (12) 

by  multiplying  both  sides  of  (11)  by  T 

= T(-Q-Un  + Sn  + 'URes(U))  (13) 

or  At 

Now  a R.unge-Kutta  method  can  be  used  to  reach  the  steady  state  in  pseudo-time.  For  stability 
reasons  the  ^ U term  is  treated  implicitly  within  the  Runge-Kuta  cycle,  for  example  the  i stage  is 
written  as, 


Qi  + a1Arr^  = «n+cvi^r(5n+7iHes'(tfi-1))  i = 1,9  (14) 

after  some  algebraic  manipulations, 
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Calculation  of  (/  + r«;Ar^;§^)  1 *n  a general  compressible  case  is  computationally  very  expen- 
sive. For  the  sake  of  simplicity,  incompressibility  is  assumed  for  the  evaluation  of  The  matrix 
equals, 
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This  simplification  seems  to  be  accurate  enough  and  has  not  brought  any  convergence  problem 
for  the  test  cases  considered  so  far. 


3 Results 

The  efficiency  of  the  dual  time-stepping  method  depends  on  the  number  of  inner  iterations  (rii ) and 
the  ratio  of  the  physical  time  step  to  the  time  step  of  a purely  explicit  scheme  (na  = ^y).  The  gain 
achieved  will  be  n2/ni-  As  a result,  the  main  goal  is  to  choose  the  physical  time  step  (At)  as  large 
as  possible  and  make  the  number  of  inner  iterations  (nj)  as  small  as  possible. 

This  section  is  divided  into  three  parts,  the  first  part  is  about  the  method  of  discretizing  the 
physical  time  derivative  and  its  effect  on  the  maximum  allowable  time  step  (At),  in  the  second  part 
the  beneficial  effect  of  preconditioning  and  multigrid  on  the  requested  number  of  inner  iterations 
is  demonstrated  and  finally  some  LES  results  of  the  channel,  cavity  and  circular  cylinder  flows  are 
presented. 
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3.1  Relation  Between  the  Time  Integration  and  the  Physical  Time  Step 

It  is  a well  known  fact  that  the  backward  differencing  scheme  unlike  the  trapezoidal  scheme  has  a 
dissipation  error  proportional  to  the  time  step.  In  order  to  see  this  effect  and  its  influence  on  the  LES 
results  a channel  flow  at  ReT  = 180  with  a mesh  of  33  x 33  x 65  points  in  the  x , y and  z directions 
is  considered.  The  time  derivative  is  discretized  either  with  a backward  differencing  (0\  = 1.5, 
00  = -2,  0- i = 0.5,  0-2  = 0.0,  7i  = 1.  and  72  = 0.0)  or  a trapezoidal  scheme  (0\  = 1.,  0o  = 0.0, 
0-1  = —1,  0-2  = 0.0,  71  = 0.5  and  72  = 0.5)  Figures  (1)  and  (2)  show  the  mean  velocity  profile 
and  the  turbulence  intensities  calculated  with  the  backward  differencing  scheme  and  the  trapezoidal 
scheme,  respectively. 

The  trapezoidal  scheme  is  much  more  flexible  than  the  backward  differencing  for  large  time  steps. 
By  increasing  the  time  step,  the  results  start  deviating  from  the  reference  solution  (obtained  with 
an  explicit  time  accurate  Runge-Kutta  method)  at  At  = 0.005  for  the  backward  differencing  and  at 
At  = 0.02  for  the  trapezoidal  scheme.  This  behavior  is  in  accordance  with  the  flat  plate  boundary 
layer  results  of  Weber  et  a/.[ll].  The  results  are  also  compared  with  the  DNS  data  of  Kim  et  a/. [12], 
the  discrepancy  between  the  LES  and  DNS  results  is  due  to  the  use  of  a second-order  scheme  used  for 
the  spatial  discretization  on  a relatively  coarse  mesh.  By  using  higher  order  methods  more  accurate 
results  cau  be  obtained,  [10]. 

3.2  The  Effect  of  Preconditioning  on  Multigrid 

The  success  of  multigrid  depends  on  the  ability  of  the  smoother  to  remove  the  high  frequency 
errors  which  cannot  be  supported  on  coarser  meshes.  For  low  Mach  number  flows  preconditioning  is 
necessary  to  remove  the  stiffness  of  the  equations  and  to  increase  the  high  frequency  damping  ability 
of  the  smoother. 

A test  was  carried  out  on  the  channel  flow  at  ReT  = 180,  M=0.06  with  a mesh  of  33  x 33  x 33 
points.  Figure  (3)  shows  part  of  the  convergence  graph  for  four  combinations  of  multigrid  and 
preconditioning.  The  number  of  inner  iterations  is  fixed  to  50. 

When  preconditioning  is  not  used  there  is  almost  no  difference  between  single  grid  and  multigrid, 
but  when  preconditioning  is  activated  multigrid  becomes  quite  efficient. 

When  both  preconditioning  and  multigrid  axe  switched  off,  the  residual  drop  of  the  first  mo- 
mentum equation  is  around  one  order  of  magnitude,  and  for  the  continuity  equation  is  less  than 
two.  When  both  preconditioning  and  multigrid  are  used,  the  residual  drops  are  more  than  four  and 
three  orders  of  magnitude  for  the  first  momentum  and  continuity  equations,  respectively.  Due  to 
the  addition  of  artificial  dissipation,  the  multigrid  works  more  efficiently  for  the  continuity  equation 
than  for  the  first  momentum  equation  in  which  the  artificial  dissipation  is  switched  off.  For  these 
calculations  a — -1  and  0 is  globally  defined  as  02  = 3.V2enl,  where  Vcmt  is  the  velocity  at  the 
centerline  of  the  channel.  A three-level  V-sawtooth  multigrid  is  used. 

3.3  Some  LES  Results 

In  all  calculations  a second-order  central  discretization  is  used  in  space  and  a second-order  trapezoidal 
scheme  in  time.  A five-stage  Runge-Kutta  scheme,  a = 0 = {1,0,0.56,0,0.44}  with 

three  evaluations  of  dissipation,  [6],  is  used  as  the  smoother  for  the  multigrid.  A three-level  V- 
sawtooth  multigrid  with  a first  order  prolongation  and  quadratic  restriction  is  used. 

3.3.1  Channel  Flow 

The  rotating  and  non-rotating  channel  flows  are  considered.  For  the  non-rotating  channel  flow,  a 
mesh  of  33  x 33  x 65  points  is  used  in  the  x,  y and  z directions,  the  streamwise,  normal  and  spanwise 
dimensions  are  47T  x 2 x 2rr.  Uniform  meshes  with  spacing  Ax~  — ~ 71  and  Az+  = ~ 18 

are  used  in  the  streamwise  and  spanwise  directions.  A non  uniform  mesh  with  cosinus  distribution  is 
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used  in  the  wall-normal  direction.  The  first  mesh  point  away  from  the  wall  is  at  y+  = ^0-  ~ 0.87 
and  the  maximum  spacing  (at  the  centerline  of  the  channel)  is  18  wall  units.  The  Reynolds  number 
based  on  the  friction  velocity  is  ReT  = 180  and  the  Mach  number  at  the  centerline  is  M = 0.06.  A 
dynamic  procedure,  [7],  is  used  to  calculate  the  Smagorinsky  coefficient  and  the  turbulent  Prandtl 
number.  For  this  calculation  a = — 1 and  02  = 3.U2cnter. 

For  the  rotating  channel  flow,  the  same  mesh  is  used  and  the  results  of  the  non-rotating  channel 
were  used  as  the  initial  solution.  The  rotation  number  is  equal  to  Ro  = 0.01.  The  rotation  number 
based  on  the  bulk  mean  velocity  (Um),  angular  velocity  (Q)  and  channel  half  width  ( h ) is, 


Ro  — 


2fi  /( 
Dm 


(17) 


The  turbulence  intensities  are  shown  in  Fig.  (4),  the  results  are  compared  with  the  DNS  data  of 
Kim  et  al. [12]  and  Kristoffersen  et  al.  [13]. 

The  physical  time  step  is  around  200  times  bigger  than  the  time  step  of  a purely  explicit  method, 
the  residual  drop  of  the  continuity  equation  is  fixed  to  -2.5  and  this  requires  almost  50  inner  iterations. 
The  gain  achieved  is  roughly  equal  to  ^ = 4. 


3.3.2  Cavity  Flow 

A lid  driven  cavity  flow  is  considered  as  another  example  to  test  the  ability  of  the  present  approach. 
The  flow  is  driven  by  the  top  wall,  z — h,  in  the  x direction  moving  with  a velocity  U — 1.  The 
Reynolds  number  is  Re  = '-£  = 10000  where  h is  the  dimension  of  the  cube.  A 33s  mesh  is  used  with 
a cosinus  distribution  of  the  mesh  points.  Starting  from  a stagnant  flow,  the  flow  is  fully  developed 
after  32/(/l/,  after  which  statistics  were  collected  during  a period  of  32/i/U.  The  mean  velocity 
profiles  along  the  vertical  and  horizontal  symmetry  lines  and  the  statistics  of  the  fluctuating  field 
are  shown  in  Fig.  5.  The  agreement  between  the  experiment,  [14],  and  simulation  is  satisfactory. 
A Smagorinsky  model  with  damping  near  the  wall  is  used  for  subgrid  scale  modeling.  For  this 
calculation  n = -1  and  01  = 10. U2.  The  physical  time  step  (At)  is  around  1000  times  bigger  than 
the  maximum  time  step  allowed  by  the  stability  condition  of  an  explicit  scheme  in  a compressible 
flow.  In  almost  130  inner  iterations  the  residual  of  the  continuity  equation  drops  four  orders  of 
magnitude.  As  a rough  estimation,  and  without  taking  into  account  the  extra  work  of  the  multigrid, 
the  implicit  approach  is  around  1000/130  ~ 7.5  times  faster  than  a purely  explicit  approach. 


3.3.3  Circular  Cylinder 

The  cylinder  has  a diameter  of  D and  a spanwise  length  of  7 tD  at  a Reynolds  number  of  Reo  = 3900. 
In  the  plane  normal  to  the  cylinder  axis,  an  O-type  grid  with  49  x 81  points  is  used  with  81  points 
on  the  surface  and  49  points  in  the  radial  direction,  33  grid  points  is  used  over  the  spanwise  length 
of  nD.  The  mesh  is  clustered  near  the  wall  to  ensure  y+  < 1 for  the  first  grid  point.  A Smagorinsky 
model  with  damping  near  the  wall  is  used  for  subgrid-scale  modeling.  For  this  calculation  a = — 1 
and  02  = 3 .U^..  The  flow  statistics  are  compared  with  the  experimental  data  of  Lourenco  and  Shih, 
taken  from,  [16].  Statistics  were  accumulated  over  T = OOD/f/^.  The  distribution  of  the  pressure 
coefficient  Cp  and  the  friction  coefficient  cj  on  the  cylinder  are  plotted  in  Fig.  6.  In  the  same  figure 
the  mean  streamwise  velocity  component  along  the  centerline  is  shown. 

Figure  7 shows  the  velocity  fluctuations  at  X = = —0.0806,  with  xi  the  location  where  the 

time  averaged  streamwise  velocity  component  is  zero. 

The  residual  drop  of  the  continuity  equation  is  fixed  to  -2.  which  is  obtained  after  100  inner 
iterations.  The  physical  time  step  of  the  implicit  method  is  around  400  times  bigger  than  the  time 
step  of  the  explicit  one.  The  implicit  method  is  y§§  = 4 times  faster  than  the  explicit  method. 
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4 Conclusions 

Preconditioning  and  multigrid  have  been  used  for  large-eddy  simulation  of  low  Mach  number  flows. 
A dual  time-stepping  approach  was  used  to  keep  the  time  accuracy  of  the  simulation.  The  channel 
and  cavity  flows  as  well  as  the  flow  around  a circular  cylinder  have  been  considered  to  test  the  ability 
of  the  method.  It  was  shown  that  the  present  method  was  4 to  7 times  faster  than  a purely  explicit 
approach.  The  efficiency  varies  according  to  the  test  case,  grid  clustering,  physical  time  step  and 
requested  residual  drop. 
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Figure  1:  Effect  of  the  time  step  with  the  backward-differencing 


Figure  2:  Effect  of  the  time  step  with  the  trapezoidal  scheme 


Figure  3:  The  effect  of  preconditioning  on  multigrid  for  the  channel  flow  calculation  with  a fixed 
number  of  inner  iterations,  left:  continuity  equation,  right:  first  momentum  equation 
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Figure  4:  Turbulence  intensities,  left:  non  rotating  channel  flow,  right:  rotating  channel  flow 


Figure  5:  (left):  mean  velocity  profile  along  two  symmetry  lines,  (middle):  lOv/uV/C/,  (right): 
10\/u/uj' /U , lines  are  used  for  the  simulation  data  and  markers  for  experimental  data. 
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Figure  6:  (left):  Cp  distribution,  (middle):  cj  distribution,  o LES  of  Breuer  [17]  (right):  mean 
streamwise  velocity,  Cp  and  < u > are  compared  with  the  experiment  of  Lourenco  and  Shih  [1G] 
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